Efficient algorithm to study interconnected networks 
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Interconnected networks have been shown to be much more vulnerable to random and targeted 
failures than isolated ones, raising several interesting questions regarding the identification and 
mitigation of their risk. The paradigm to address these questions is the percolation model, where 
the resilience of the system is quantified by the dependence of the size of the largest cluster on 
the number of failures. Numerically, the major challenge is the identification of this cluster and 
the calculation of its size. Here, we propose an efficient algorithm to tackle this problem. We 
show that the algorithm scales as 0(N log N), where N is the number of nodes in the network, a 
significant improvement compared to 0(N 2 ) for a greedy algorithm, what permits studying much 
larger networks. Our new strategy can be applied to any network topology and distribution of 
interdependencies, as well as any sequence of failures. 
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I. INTRODUCTION 

Most real networks are strongly dependent on the func- 
tionality of other networks [TH5]. For example, the per- 
formance of a power grid is assured by a system of global 
monitoring and control, which depends on a communica- 
tion network. In turn, the servers of the communication 
network rely on the power grid for power supply. This 
interdependence between networks strongly affects their 
resilience to failures. Buldyrev et al. [%\ have developed 
the first strategy to analyze this coupling in the frame- 
work of percolation. To the conventional representation 
of complex networks, where nodes are the agents (e.g., 
power stations or servers) and edges are the interactions 
(either physical or virtual), they added a new type of 
edges, namely, dependency links, to represent the inter- 
network coupling. Such links couple two nodes from dif- 
ferent networks in such a way that if one fails the other 
cannot function either. They have shown that this cou- 
pling promotes cascading failures and strongly affects the 
systemic risk, drawing the attention towards the dynam- 
ics of coupled systems. A different framework based on 
epidemic spreading has also been proposed leading to the 
same conclusions [7J. 

To quantify the resilience of interconnected networks, 
one typically simulates a sequence of node failures (by re- 
moving nodes) and measures the dependence of the size 
of the largest connected component on the number of fail- 
ures 0[S]- The first studies have shown that, depending 
on the strength of the coupling (e.g., fraction of depen- 
dency links), at the percolation threshold, this function 
can change either smoothly (weak coupling) or abruptly 
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(strong coupling) [TT]. As reviewed in Refs. [TUHT^] . sev- 
eral works have followed studying, for example, the de- 
pendence on the coupling strength [51 [IS] , the role of net- 
work topology, and the phenomenon on geographically 
embedded networks [141 115] . A more general framework 
was also developed to consider a network of networks 
[TrjHTB] . In all cases, astonishing properties have been re- 
vealed, which were never observed for isolated systems. 

For many cases of interest, the size of the largest com- 
ponent needs to be computed numerically as the available 
analytic formalisms are limited to very simple networks, 
interdependencies, and sequence of failures EJ fT3] , 
However, the determination of this largest component 
and its size is not a trivial task. When a node is removed 
(fails), the triggering of cascading failures and multiple 
interdependencies need to be considered. Here we pro- 
pose an efficient algorithm, where a special data structure 
is used for the fast identification of the largest fragment 
when the network breaks into pieces. We show that the 
algorithm scales as 0(N log N) , where N is the number of 
nodes in the network, while the one of a greedy algorithm 
is 0{N 2 ). This strategy permits studying very large sys- 
tem sizes and many samples, which leads to much more 
accurate statistics. Since our description is generic, it is 
possible to consider any network and distribution of in- 
terdependencies, as well as sequences of failures [HH5T]. 

The paper is organized in the following way. The algo- 
rithm is described in Sec. [H] and its efficiency discussed 
in Sec. Mil In Sec. IIVI we make some final remarks and 
discuss possible future applications. 



II. ALGORITHM 

Figure [l] shows the dependence of the fraction of nodes 
in the largest connected cluster s on the fraction of re- 
moved nodes 1 — p, for two Erdos-Renyi networks with 



2 





FIG. 1. (Color online) Comparison between numerical and 
analytic results for percolation on two coupled Erdos-Renyi 
networks. In the main plot, the numerical (data points) re- 
sults for the size of the largest connected cluster are obtained 
from two coupled networks with N = 512000 nodes each and 
coupling strengths of q = 1,0.9, and 0.7 (from left to right). 
The lines correspond to the analytic results computed as in 
Ref. [22]. The numerical results are averages over 50 differ- 
ent pairs of networks and 100 sequences of random failures. 
All numerical and corresponding analytic curves overlap. In 
the inset we see the number of iterations per failure, defined 
as the number of times a cascade triggered by a node failure 
propagates back and forth between the two networks. The 
lines tag the analytic results for the percolation threshold. 
Initially, a node removal does not trigger any cascade but, as 
one approaches the percolation threshold, very large cascades 
occur, resulting in the collapse of the entire system. 



more than one million nodes. When nodes are randomly 
removed, the largest connected component decreases in 
size, until the network is completely fragmented above a 
threshold 1 — p c . In the inset, we see the evolution of 
the number of iterations per removed node. An iteration 
corresponds to a set of failures in one network triggered 
by an internetwork coupling, i.e., by the removal of a 
dependency link. The number of iterations is negligibly 
small for low values of 1 — p but peaks at the thresh- 
old. Following the cascade after removing a node is the 
most computational demanding task. Consequently, an 
efficient algorithm is required to identify, in a fast way, if 
a node removal triggers a cascade or not. 

Here, we propose an efficient data structure to recog- 
nize the beginning of a cascade and identify the different 
fragments resulting from a node removal. Since we are 
interested in the evolution of the largest connected com- 
ponent, we only follow this cluster. Our algorithm uses 
a hierarchical data structure with different levels. As il- 
lustrated in Fig. |2j we choose the node with the highest 
degree as the root and assign to it the level L = 0. All 
neighbors of this root are on the second level (L = 1) and 
they are directly connected to the root. All neighbors of 
the second level, which have not an assigned level yet, are 



FIG. 2. Example of the level structure for a small net- 
work. An energy number is assigned to each node, which 
corresponds to the shortest distance to the root node (the 
most central node in this example). 



then placed on the third level. We proceed iteratively in 
the same way, until all nodes of the cluster have a level. 
Note that we can have links within the same level and 
between levels but, in the latter, the level difference is 
limited to unity. The depth of the level structure is the 
maximal distance between the root and any other node 
in the network. For random networks, this depth ap- 
proximately scales with log N [21] and it scales even 
slower for many scale-free networks [2.5]. Note that, in 
the case of n coupled networks we will have n different hi- 
erarchical structures, i.e., one per network, representing 
its largest component. 

When a node in level L is removed, the ordering needs 
to be updated. All neighbors at a higher level L+l which 
are connected to another node in level L remain in the 
same level, as shown in Fig. [3j The nodes in level L + l 
which have no further neighbors in level L but only in 
level L + l, need to be updated (moved one level up) as 
well as the entire branch connected to them. In those two 
cases, the size of the largest connected component in this 
iteration is just changed by unity (the initially removed 
node). If neither of those cases occurs, i.e. all neighbors 
have a higher level, we proceed iteratively through the 
branch of neighbors with a breadth first search (up in 
level) until we detect one node in level L' which has at 
least one neighbor in level L' or L' — 1 which is not de- 
tected by the breadth first search. In this case, the entire 
branch of detected nodes is updated, starting from the 
last node in level V . On the other hand, if no node in the 
branch establishes a connection with the other branches, 
it implies that the largest component was split into sub- 
networks and one has to decide which one is the largest. 
Then the size of the largest connected component is ad- 
justed and all nodes reorganized (see example in Fig. [4j) . 



III. NUMBER OF COMMANDS AND 
COMPUTATIONAL TIME 

To assess the efficiency of the algorithm, we study the 
dependence of the number of commands Cat on the net- 
work size N. We count as a command, every time a node 
in one of the networks is removed, its level changed, or 
just checked during the reorganization of the level struc- 
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FIG. 3. (Color online) Example of a node removal. In 
this case, the red (light) node with L = 1 is removed. Since 
the only neighbor of this node with higher level has another 
neighbor with L = 1, the size of the largest connected cluster 
is only reduced by one. 






I l 
I I 
I I 

o ! l ! 2 

FIG. 4. (Color online) Example of a reorganization due to a 
node removal. In the example, the root node is removed, thus 
the level structure has to be reorganized. First, the largest 
subnetwork is identified and all the other subnetworks are 
removed. One of the surviving first neighbors of the old root 
is randomly selected to become the new root. Then the levels 
of its neighbors are updated as well as their branches. Note 
that the update of a branch is complete when the level of 
a node remains the same. In the worst case scenario, the 
complexity of the entire update process is O(M), where M is 
the size of the entire branch. 



ture. Figure [5] (main plot) shows the size dependence of 
the average Cn for two coupled Erdos-Renyi networks 
[26] . The line fitting the data points is C \og(N)N Cl , 
where C\ is 1.02 ± 0.03. In a greedy algorithm where 
the largest connected cluster is recalculated by counting 
all remaining nodes in this component after each node 
removal, the number of commands is expected to scale 
as 0(N 2 ). With our data structure, this limit where all 
nodes are checked, would correspond to the worst case 
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FIG. 5. (Color online) The number of used commands in the 
program Cjv versus the system size N. The function f(N) = 
C log(N)N Cl is fitted to the observed values. The fitting 
parameter C\ is 1.02 ± 0.03 for Erdos-Renyi networks and 
0.97 ± 0.04 for scale-free networks, respectively. 



scenario, where the removed nodes would systematically 
be the root. Therefore, our algorithm represents a signif- 
icant improvement over the traditional greedy algorithm. 
Also in Fig. [5j we plot the number of commands Cm for 
two coupled scale-free networks, with degree exponent 
7 = 2.5 [57] • The same scaling with the network size was 
found, with C\ = 0.97 ±0.04. So, for the two considered 
types of networks, our algorithm scales with O(NlogN). 
In general, this scaling will depend on how the average 
shortest path scales with the network size. 

Figure [6] shows the size dependency of the average com- 
putational time t(N) required to compute an entire se- 
quence of node removals. We show the ratio t(N)/t(N/2) 
obtained from two computers with 6 MB and 12 MB 
CPU cache for Erdos-Renyi networks (main plot) and 
scale-free networks (inset). In both cases, we observe a 
crossover between two different scaling regimes at a cer- 
tain system size N* . This crossover at N* — 4000 and 
N* = 8000 for 6 MB cache and 12 MB cache, respec- 
tively, depends on the size of the CPU cache memory 
(L2). For network sizes N < N* , the size of the system 
is such that all information can be kept inside the CPU 
cache, being more efficient. For N > N* , not all informa- 
tion fits in the CPU cache and the efficiency decreases, 
since the access to the Random Access Memory (RAM) 
is slower. 



In the first regime N < N* the increase of the 
CPU time is consistent with an algorithm scaling as 
0(N log N). In the second regime N > N* the CPU 
time seems to converge to the same logarithmic scaling. 
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FIG. 6. System-size dependence of the average CPU time 
t(N), in seconds, necessary to calculate one sequence of node 
removals. We see the ratio between two times t(N) /t(N/2) for 
two different machines (different capacity of the CPU cache 
memory) for two coupled Erdos-Renyi (ER, main plot) and 
scale-free networks (SF, inset). Results are averages over 100 
sequences of random removals and 50 different networks. 



dent networks. Although, in general, the algorithm can 
be considered to study percolation in both isolated and 
coupled networks, its is tailored for coupled ones. We 
have shown that the algorithm complexity is 0(N log N) , 
a significant improvement over the greedy algorithm of 
complexity 0(N 2 ). 

With our efficient algorithm, it is now possible to sim- 
ulate much larger system sizes, a relevant feature to de- 
velop accurate studies. One of the most striking results 
of coupled networks is that, for strong coupling, the frag- 
mentation of the network into pieces occurs in a discon- 
tinuous way [6]. The possibility of accurate measure- 
ments, with reduced finite-size effects, permits to de- 
termine with high precision the critical coupling above 
which the percolation transition is discontinuous |13) . 
Our algorithm can now be applied to any network topol- 
ogy, sequence of failures (e.g., random or high-degree 
node) [5] , and distribution of dependency links [33] , help- 
ing clarifying how to mitigate the systemic risk stemming 
from interdependencies. 
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